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Abstract 

The  high  order  spectral  element  approximation  of  the  Euler  equations  is  stabilized  via  a 
dynamic  sub-grid  scale  model  (Dyn-SGS).  This  model  was  originally  designed  for  linear  finite 
elements  to  solve  compressible  flows  at  large  Mach  numbers.  We  extend  its  application  to 
high-order  spectral  elements  to  solve  the  Euler  equations  of  low  Mach  number  stratified  flows. 
The  major  justification  of  this  work  is  twofold:  stabilization  and  large  eddy  simulation  are 
achieved  via  one  scheme  only. 

Because  the  diffusion  coefficients  of  the  regularization  stresses  obtained  via  Dyn-SGS  are 
residual-based,  the  effect  of  the  artificial  diffusion  is  minimal  in  the  regions  where  the  solution 
is  smooth.  The  direct  consequence  is  that  the  nominal  convergence  rate  of  the  high-order 
solution  of  smooth  problems  is  not  degraded.  To  our  knowledge,  this  is  the  first  application  in 
atmospheric  modeling  of  a  spectral  element  model  stabilized  by  an  eddy  viscosity  scheme  that, 
by  construction,  may  fulfill  stabilization  requirements,  can  model  turbulence  via  LES,  and  is 
completely  free  of  a  user-tunable  parameter. 

From  its  derivation,  it  will  be  immediately  clear  that  Dyn-SGS  is  independent  of  the  nu¬ 
merical  method;  it  could  be  implemented  in  a  discontinuous  Galerkin,  finite  volume,  or  other 
environments  alike.  Preliminary  discontinuous  Galerkin  results  are  reported  as  well.  The 
straightforward  extension  to  non-linear  scalar  problems  is  also  described.  A  suite  of  ID,  2D, 
and  3D  test  cases  is  used  to  assess  the  method,  with  some  comparison  against  the  results 
obtained  with  the  most  known  Lilly-Smagorinsky  SGS  model. 


1  Introduction 

The  search  for  the  best  stabilization  method  for  high-order  continuous  Galerkin  (CG)  to  solve  the 
Euler  equations  remains  an  active  field  of  research  [24] .  This  is  partially  due  to  certain  properties 
that  the  stabilization  method  should  have  but  that  are  not  always  present.  In  particular,  if  achieved 
via  some  type  of  artificial  diffusion,  (i)  the  effect  of  stabilization  on  the  solution  should  vanish 
where  the  solution  is  smooth* 1  and  (ii)  should  decrease  as  the  grid  is  refined.  Most  stabilized 
(Galerkin)  methods  possess  the  latter  but  not  the  first  (e.g,  [21],  spectral  elements).  To  fulfill 
the  first  requirement,  solution-dependent  viscosities  have  existed  for  quite  some  time,  especially  to 
stabilize  low  order  finite  element  [7;  29;  28;  45;  27],  finite  volume  [1],  and  high-order  discontinuous 
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1This  property  is  referred  to  as  consistency. 


1 


Report  Documentation  Page 


Form  Approved 
OMB  No.  0704-0188 


Public  reporting  burden  for  the  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing  data  sources,  gathering  and 
maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this  collection  of  information, 
including  suggestions  for  reducing  this  burden,  to  Washington  Headquarters  Services,  Directorate  for  Information  Operations  and  Reports,  1215  Jefferson  Davis  Highway,  Suite  1204,  Arlington 
VA  22202-4302.  Respondents  should  be  aware  that  notwithstanding  any  other  provision  of  law,  no  person  shall  be  subject  to  a  penalty  for  failing  to  comply  with  a  collection  of  information  if  it 
does  not  display  a  currently  valid  OMB  control  number. 


1.  REPORT  DATE 

JAN  2015 


2.  REPORT  TYPE 


4.  TITLE  AND  SUBTITLE 

Stabilized  High-order  Galerkin  Methods  Based  on  a  Parameter-free 
Dynamic  SGS  Model  for  LES 

6.  AUTHOR(S) 


7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

Naval  Postgraduate  School, Department  of  Applied 
Mathematics, Monterey, CA, 93943 

9.  SPONSORING/MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 


3.  DATES  COVERED 

00-00-2015  to  00-00-2015 

5a.  CONTRACT  NUMBER 

5b.  GRANT  NUMBER 

5c.  PROGRAM  ELEMENT  NUMBER 

5d.  PROJECT  NUMBER 

5e.  TASK  NUMBER 

5f.  WORK  UNIT  NUMBER 

8.  PERFORMING  ORGANIZATION 
REPORT  NUMBER 


10.  SPONSOR/MONITOR'S  ACRONYM(S) 

11.  SPONSOR/MONITOR'S  REPORT 
NUMBER(S) 


12.  DISTRIBUTION/AVAILABILITY  STATEMENT 

Approved  for  public  release;  distribution  unlimited 

13.  SUPPLEMENTARY  NOTES 

Journal  of  Computational  Physics  01/2015;  301:77-101 

14.  ABSTRACT 

The  high  order  spectral  element  approximation  of  the  Euler  equations  is  stabilized  via  a  dynamic  sub-grid 
scale  model  (Dyn-SGS).  This  model  was  originally  designed  for  linear  finite  elements  to  solve  compressible 
flows  at  large  Mach  numbers.  We  extend  its  application  to  high-order  spectral  elements  to  solve  the  Euler 
equations  of  low  Mach  number  stratified  flows.  The  major  justification  of  this  work  is  twofold: 
stabilization  and  large  eddy  simulation  are  achieved  via  one  scheme  only.  Because  the  diffusion  coefficients 
of  the  regularization  stresses  obtained  via  Dyn-SGS  are  residual-based,  the  effect  of  the  artificial  diffusion 
is  minimal  in  the  regions  where  the  solution  is  smooth.  The  direct  consequence  is  that  the  nominal 
convergence  rate  of  the  high-order  solution  of  smooth  problems  is  not  degraded.  To  our  knowledge,  this  is 
the  first  application  in  atmospheric  modeling  of  a  spectral  element  model  stabilized  by  an  eddy  viscosity 
scheme  that  by  construction,  may  fulfill  stabilization  requirements,  can  model  turbulence  via  LES,  and  is 
completely  free  of  a  user- tunable  parameter.  From  its  derivation,  it  will  be  immediately  clear  that 
Dyn-SGS  is  independent  of  the  numerical  method;  it  could  be  implemented  in  a  discontinuous  Galerkin, 
finite  volume,  or  other  environments  alike.  Preliminary  discontinuous  Galerkin  results  are  reported  as 
well.  The  straightforward  extension  to  non-linear  scalar  problems  is  also  described.  A  suite  of  ID,  2D  and 
3D  test  cases  is  used  to  assess  the  method,  with  some  comparison  against  the  results  obtained  with  the  most 
known  Lilly- Smagorinsky  SGS  model. 

15.  SUBJECT  TERMS 


16.  SECURITY  CLASSIFICATION  OF: 

17.  LIMITATION  OF 

18.  NUMBER 

19a.  NAME  OF 

ABSTRACT 

OF  PAGES 

RESPONSIBLE  PERSON 

a.  REPORT 

unclassified 

b.  ABSTRACT 

unclassified 

c.  THIS  PAGE 

unclassified 

Same  as 
Report  (SAR) 

38 

Standard  Form  298  (Rev.  8-98) 

Prescribed  by  ANSI  Std  Z39-18 


Galerkin  (DG)  [44;  34]  methods.  Residual-based  viscosity  for  spectral  elements  is  found  in  the 
extensive  work  on  the  entropy- viscosity  method  by  [22;  23].  Based  on  the  idea  of  the  Spectral 
Vanishing  Viscosity  of  [50],  the  entropy- viscosity  technique  was  designed  in  such  a  way  that  (iii) 
the  spectral  accuracy  of  the  spatial  approximation  is  preserved  if  the  solution  is  smooth.  This 
represents  the  third  property  that  a  viscosity-based  stabilizing  scheme  should  possess. 

A  computationally  inexpensive  sub-grid  scale  (SGS)  model  for  compressible  Large  Eddy  Sim¬ 
ulation  (LES)  that  possesses  (i)  and  (ii)  was  recently  designed  to  stabilize  linear  finite  elements 
for  different  flow  regimes  [43].  At  first  sight,  this  model  is  a  close  relative  of  the  entropy- viscosity 
method,  with  the  fundamental  difference  that  no  entropy  equation  is  used  to  construct  the  dynamic 
dissipation.  Rather,  the  residuals  of  the  Euler  equations  are  used  to  build  fully  coupled  dissipation 
coefficients.  Furthermore,  unlike  the  entropy-viscosity  method  and,  to  the  best  of  our  knowledge, 
unlike  any  of  the  dissipation  models  listed  above,  this  scheme  converges  to  the  unique  entropy 
solution  (see  the  proof  in  [42]).  The  use  of  the  residuals  of  each  equation  in  the  system  makes 
the  artificial  diffusion  non-linear.  The  reason  for  exploring  non-linear  dissipation  stems  from  the 
analysis  reported  in  [15],  where  it  is  concluded  that  linear  stabilization  can,  at  most,  give  a  solution 
that  converges  to  a  weak  solution  that  is  not  the  entropy  solution. 

In  the  current  paper,  we  explore  the  capabilities  of  this  dynamic  model  (from  now  on,  Dyn-SGS) 
for  the  stabilization  of  the  spectral  element  solution  of  the  Euler  equations  for  low  Mach  number 
stratified  flows.  Furthermore,  we  extend  the  method  for  the  solution  of  quasi-linear  scalar  equations 
for  use  in  transport  schemes  as  required  in  atmospheric  modeling. 

The  derivation  of  Dyn-SGS  is  tied  to  the  fundamental  idea  of  LES;  the  flow  equations  (Euler, 
in  this  particular  case)  are  filtered  to  separate  the  large  (grid  resolved)  from  the  small  (sub-grid, 
unresolved)  scales  [35].  The  filtering  operation  yields  a  new  set  of  equations  containing  additional 
terms  that  account  for  the  effect  of  the  unresolved  scales  onto  the  resolved  solution. 

The  original  version  of  the  current  method  [43]  relies  on  the  regularization  of  the  continuity 
equation,  other  than  momentum  and  energy.  The  idea  of  mass  stabilization  is  not  new;  we  find  it  in, 
e.g.,  [44;  34;  49;  5;  23;  41;  24].  Despite  the  remarkable  results  that  are  shown  in  the  aforementioned 
literature,  this  type  of  regularization  is  often  the  subject  of  criticism  by  physicists  who,  for  the  most 
part,  doubt  the  physical  meaning  of  mass  dissipation  on  the  one  hand  and  the  conservation  of  mass 
on  the  other.  To  not  affect  conservation,  in  the  current  work  we  build  Dyn-SGS  so  that  dissipation 
is  not  applied  to  the  continuity  equation.  Its  formal  derivation  is  described  in  the  text.  Note  that 
mass  may  still  be  conserved  when  the  continuity  equation  is  stabilized,  as  long  as  certain  boundary 
conditions  are  applied. 

In  Dyn-SGS  the  diffusion  coefficients  of  the  (turbulent)  stresses  are  a  function  of  the  equation 
residuals.  In  contrast  to  most  artificial  diffusion  methods  that  apply  the  same  operator  and  coeffi¬ 
cient  to  every  equation,  Dyn-SGS  applies  separate  diffusion  coefficients  to  each  equation.  Because 
Dyn-SGS  is  residual-based,  the  effect  of  the  artificial  diffusion  is  minimal  in  the  regions  where  the 
solution  is  smooth.  The  third  property  (iii)  mentioned  above  is  a  direct  consequence  of  this. 

In  spite  of,  e.g.,  its  higher  parallel  cost  than  its  second  order  counterpart  and  the  possible  need  of 
a  tunable  parameter,  linear  and  monolithic2  hyper-viscosity  remains  the  classical  approach  used  by 
atmospheric  modelers  [31].  A  recent  review  of  dissipation-based  stabilization  schemes  for  Galerkin 
methods  can  be  found  in  [40].  This  work  represents,  to  our  knowledge,  the  first  application  of  a 
non-linear  and  grid  dependent  artificial  diffusion  to  stabilize  high-order  spectral  elements  in  the 
solution  of  the  Euler  equations  for  stratified  flows  and  that  possesses  the  three  properties  (i),  (ii), 

2 Monolithic  [24]  means  that  the  2o-order  operator  (— 1)“+1V“  •  (/iV“)  with  diffusion  coefficient  //  is  applied 
equally  to  all  quantities. 
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Table  1:  List  of  constants  used  in  this  paper 


Symbol 

Description 

Value 

cp 

specific  heat  of  dry  air  at  constant  pressure 

1004.5  Jkg^1  K— 1 

cv 

specific  heat  of  dry  air  at  constant  volume 

717.5  Jkg-1  K_1 

7 

Cp/  cv 

1.4 

R 

gas  constant  of  dry  air 

287.17  Jkg-1  K— 1 

9 

magnitude  of  acceleration  of  gravity 

9.80616  ms-1 

po 

reference  pressure 

105  Pa 

uj 

angular  velocity  of  the  Earth 

7.292  x  10-5  s_1 

and  (iii).  Moreover,  Dyn-SGS  is  fully  parameter  free.  The  extension  to  linear  and  non-linear  scalar 
problems  is  straightforward  and  is  described  as  well. 

The  current  method  is  successfully  tested  on  a  series  of  ID,  2D  and  3D  benchmark  problems 
for  stratified  and  gravity  driven  atmospheres  and  scalar  transport  problems. 

This  study  is  a  first  effort  towards  the  construction  of  an  LES  model  within  the  Nonhydrostatic 
Unified  Model  of  the  Atmosphere  (NUMA)  developed  by  the  authors  during  the  past  few  years 
[33;  19]. 

The  rest  of  the  paper  is  organized  as  follows.  The  sets  of  equations  and  the  SGS  model  are 
described  in  Section  2.  The  numerical  discretization  method  of  these  equations  is  reported  in  Section 
3.  The  results  for  the  Euler  and  scalar  equations  are  reported  in  Sections  4  and,  respectively,  5. 
The  conclusions  are  given  in  Section  6. 


2  Governing  equations 


Let  SlgR3  be  a  fixed  3-dimensional  domain  with  boundary  T  and  Cartesian  coordinates  x  = 
[xi(x),X2(y),X3(z)].  Let  us  identify  the  dry  air  density,  the  velocity  vector,  and  the  potential 
temperature  with  the  symbols  p1  u  =  iq  (i  =  1, 2, 3)  and  6.  Then,  the  time-dependent  Euler  equations 
in  advective  form  can  be  written  as: 


dp  dpuj_  = 
dt  dxj 

dui  dui  1  dp 
dt  +  tjl  dxj  +  p  dxi 


-gdi, 


(la) 

(lb) 


d9 

dt 


dO 

'  dxj 


=  0, 


(lc) 


where  g  is  the  acceleration  due  to  gravity  that  acts  in  the  direction  of  £1,2,3  =[0  0  —  1].  Equations 
(1)  must  be  solved  in  0  Vt  £  R-1"  given  proper  initial  and  boundary  conditions.  Pressure,  p,  is 
related  to  0 ,  and  p  through  the  equation  of  state  for  a  perfect  gas 


P  =  P  0 


The  values  of  the  reference  pressure,  po ,  and  of  the  other  constants  used  in  the  paper  are  reported 
in  Table  1. 
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Remark  1  Because  in  Numerical  Weather  Prediction  (NWP)  the  small  variations  of  the  ther¬ 
modynamic  quantities  ( p,0,p )  with  respect  to  a  constant  and  hydrostatically  balanced  reference 
background  state  (pref^ref^Pref)  are  the  quantities  of  major  interest  [38],  the  splitting  A p(t,x)  = 
p(f,x)  —  pref(z)1  A0(f,x)  =  0(f,x)  —  9ref{z ),  and  A p(f,x)  =  p(f,x)  —  pref(z)  must  be  introduced  to 
solve  Eqs.  (1).  Nevertheless,  to  keep  the  notation  clear,  the  symbol  A  will  only  be  used  in  the 
analysis  of  the  results  rather  than  in  the  definition  and  analysis  of  the  LES  equations  described 
from  now  on. 

2.1  LES  and  the  Dyn-SGS  model 

The  LES  formulation  is  obtained  by  first  introducing  the  spatial  filtering  operation 

7(x)  =  [  G-K(x-x)f{x)dX,  (2) 

Jn 

where  the  resolved  quantity,  /,  is  obtained  from  the  filtering  function  G  of  the  instantaneous 
quantities,  /,  using  a  filter  width  A.  The  application  of  (2)  to  the  continuity  equation  (la)  results 
in  the  presence  of  an  additional  sub-grid  term  on  the  right-hand  side.  To  avoid  it,  the  change  of 
variable  f  =  pf  /p  [16]  is  also  introduced3.  The  details  of  the  filtering  operations  for  the  compressible 
equations  can  be  found  in,  e.g.  [18]  and  citations  therein.  The  two  operations  yield  the  filtered 
equations 


dp  dpuj  _ 
dt  dxj 

duj  ~  duj  1  dp  _  1  dr^GS 
dt  J  dxj  p  dxi  p  dxj  9 

dO  _  dO  _  1  dQS°s 

dt  dxj  p  dxj  ’ 


(3a) 

(3b) 

(3c) 


where  the  derivatives  of  rfjGS  and  Q^GS  on  the  right-hand  side  of  (3b)  and  (3c)  represent  the 
contribution  of  the  unresolved  sub-grid  scales  (SGS).  In  (3b),  t^gs  is  the  turbulent  stress  tensor, 


T-j  =  p{UiUj  -UiUj)  , 


which  is  modeled  as  a  function  of  the  tensor 

_  1  f  duj  duj\ 

lJ  2  \dxj  dxi )  ’ 

as 

T?GS=2pSij.  (4) 

How  the  coefficient  p  is  constructed  defines  the  SGS  method  at  hand.  It  will  be  defined  shortly. 
Similarly,  the  quantity  QjGS  in  (3c)  results  from  filtering  Eq.  (lc);  it  is  given  by 

Qj°S  =  P  (®uj  ~  (5) 

3 From  now  on,  the  symbols  7  and  •  will  indicate  the  grid-resolved  quantities  (large-eddy). 
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Table  2:  Element  characteristic  length  for  anisotropic  grids 


Formula  Notes 

A  =  (AxAyAz)1/3  /p 
A  =  f  (AxAyAz)1^  /p 

where  /  =  cosh  (in2  ai  +  In2  a2  —  lnai  In  <22)]  , 

where  a\ ,2  —  A smalli  2/ max  (A#,  Ay,  Az) /(AT  +  1)  A Smalli  2  are  the  two  smallest  sizes. 

N  is  the  order  of  approximation.  See  Section  3. 


Reference 

[13] 

[46] 


and  is  modeled  as 


QfGS 


(6) 


where  the  coefficient  k  will  be  defined  below.  To  proceed  with  the  construction  of  p  and  k,  we 
introduce  the  residuals  of  the  filtered  equations  from  which  the  dissipative  operators  are  excluded. 
We  write: 


R(p) 


dp  dpuj 
dt  dxj 


(7a) 


R(ui) 


dui  _  dui  1  dp 
dt  J  dxj  p  dxi 


gSi, 


(7b) 


n,~  de  _  do 

m~~dt+Ujdx 7 


(7c) 


The  time  derivative  in  (7)  can  be  approximated  via  a  finite  difference  approximation.  As  previously 
done  in,  e.g.,  [23]  we  use  a  2nd-order  backward  differentiation  formula.  Other  methods  can  be  used. 
The  coefficients  p  and  n  are  calculated  element-wise  on  every  high  order  element  Oe  given  a  spectral 
element  approximation  of  Equations  (3)  -  the  spectral  element  method  will  be  described  in  Section 
3  -  In  the  current  model,  the  filter  width  is  taken  as  the  characteristic  size  of  an  element.  Simply, 
given  an  element  f le  of  order  N  and  edge  lengths  Ax,Ay,Az  of  comparable  orders  of  magnitude, 
we  define  the  following  characteristic  length,  and  hence  filter  width: 


A  =  min  (Ax,Ay,Az)  /  (N  + 1). 


This  definition  is  sufficient  given  the  scope  of  the  current  study;  nevertheless,  a  more  proper  defi¬ 
nition  of  A  for  LES  should  be  used  in  future  work. 


Remark  2:  selection  of  the  element  characteristic  length  The  choice  of  the  characteristic 
element  size  is  tied  to  the  aspect  ratio  of  the  element  at  hand.  For  grids  made  of  relatively  reg¬ 
ular  hexahedra,  being  Ax,  Ay  and  Az  of  the  same  order  of  magnitude,  the  choice  is  trivial.  On 
the  contrary,  highly  anisotropic  grids  (usually  the  case  for  global  circulation  models  of  the  Earth 
atmosphere,  or  boundary  layer  grids),  stabilization  may  be  greatly  affected  by  an  improper  choice 
of  A.  Based  on  our  experience,  for  Az/(Ax,Ay)  <C  1  and  Ax/ Ay  ss  1,  we  recommend  to  split  the 
dynamic  diffusion  operator  into  a  horizontal  and  a  vertical  component.  In  this  way,  simulations 
on  anisotropic  grids  will  rely  on  anisotropic  diffusion.  Should  a  splitting  operation  not  be  possible, 
then  we  recommend  the  options  in  Table  2. 
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For  the  sensible  temperature  T  =  0(p/po)R^Cp  and  one  element  of  characteristic  length  A,  we 
start  by  defining  the  following  quantities: 


Mmax|ne  —  0.5A  — 


|S|  + 


(8) 


where  |u|  +  \fjpjp  is  the  maximum  wave  speed  and 


M res|ne  =  A"  max 


(  ll^(d)l|oo,ne 

\  lld-p||oo,n 


||-R(u)||oo,ne  H-RWIIoo.qA 
||S  — u||oo,n  ’  He-eilocn  ) 


(9) 


where  ^  indicates  the  space  average  of  the  quantity  at  hand  over  ft  and  the  norms  ||  •  ||oo,n  at  the 
denominator  are  used  for  normalization  to  preserve  the  correct  dimension  of  the  resulting  equation. 
Having  pmax  and  pres  constructed,  we  can  compute  the  dynamic  coefficients  of  the  viscosity  terms 
as: 

A*  lne  =  min  (/rmax|ne ,  A*res  |oe)  (10a) 

«|f2e  =  (10b) 

where  Pr  is  an  artificial  Prandtl  number  whose  value ^will  be  defined  in  Section  4. 

The  size  of  the  residuals  is  proportional  to  (Aqe)  at  discontinuities  and  is  relatively  small  near 
smooth  regions.  Moreover,  the  size  of  the  viscosity  terms  never  exceeds  /imax)  which  is  equivalent 
to  a  stable  upwind  scheme  in  regions  with  sharp  discontinuities.  We  briefly  illustrate  this  last  point 
by  taking  a  finite  difference  approximation  of  the  ID  advection-diffusion  equation  with  positive 
velocity  u  and  a  diffusion  coefficient  proportional  to  pmax-  Given  n  time  steps  and  a  grid  of  j 
points,  we  write  the  following: 


dq  dq  1  —  d2q 

5i+“ai~2a“aS5“ 

A  t  2A  2  A^2 

At  +2U  A 

n+i  n  n_  n 

huQj 


At  A 

where  the  last  expression  is  an  upwinded  difference  approximation. 


Remark  3:  a  parallel  with  Lilly- Smagorinsky  Going  back  to  Eq.  (4)  and  the  definition  of 
/i,  the  eddy  viscosity  model  due  to  Lilly  and  Smagorinsky  [37;  47]  (valid  for  low  Mach  number  flows 
only)  reads  as  follows: 

M  =  2C'sA2|S|v/r:rRi,  (12) 
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where  Cs  is  a  constant  to  be  determined  (see  [37;  36]),  |S|  =  ^/2 SijSij,  and  Ri  is  the  Richardson 
number  defined  as 


Ri 


g_M 
6  dz 


/|S|- 


Given  /i,  k  is  obtained  analogously. 


Dyn-SGS  and  Smagorinsky  clearly  come  from  the  same  root.  This  means  that  Dyn-SGS  can  be 
implemented  with  ease  within  any  code  that  already  contains  an  implementation  of  Smagorinsky. 
By  theoretical  means  and  physical  assumptions,  the  Smagorinsky  constant  assumes  values  Cs  = 
0.2  —  0.22;  these  are  twice  as  large  as  the  values  needed  for  LES  calculations  of  practical  importance 
[9].  The  more  typical  value  Cs  =  0.14  is  used  in  the  Smagorinsky  runs  reported  in  Section  4.  No 
user-defined  parameter  has  to  be  defined  for  Dyn-SGS. 


2.2  Extension  to  quasi-linear  scalar  equations 

The  extension  to  scalar  problems  is  straightforward.  We  illustrate  it  for  the  three-dimensional 
advection  equation  of  a  scalar  q: 

<*> 


where  Uj  are  the  components  of  a  prescribed  velocity  field.  When  q  =  u\,  and  j  =  1,  the  ID  quasi- 
linear,  inviscid  Burgers  equation4  [4;  8]  is  recovered.  If  we  filter  (13)  as  done  for  the  Euler  equations, 
we  find  the  new  problem 


dq  d  q  dQ^GS 

dt  dxj  dxj 


(14) 


As  in  Eq.  (6),  Q^GS  is  modeled  as 


QfGS 


where  fi  is  the  residual-dependent  function 


s\ne 


0.5A2 


P(g)l|oo,ne 

lk-9lloo,n  ’ 


(15) 


3  Variational  formulation  and  spectral  element  approxima¬ 
tion 

Equations  (3)  are  solved  via  the  spectral  element  method  in  the  domain  12  with  boundary  T.  To 
proceed,  we  define  the  following  notation:  given  the  Sobolev  space  Hl{  12)  of  weakly-differentiable 
functions,  the  space  V  C  H1  of  test  and  trial  functions  of  the  Galerkin  formulation  is  introduced. 
Given  the  space  L2  of  real- valued  functions  that  are  square  integrable  in  12,  the  2-norm  associated 
with  it  is  denoted  by  ||  •  || 2-  Given  a  finite  element  partition  Qh  =  U”f.112e  of  the  domain  12  into  ne 
conforming  hexahedra  12e  of  characteristic  length  A,  Vh  is  the  finite  dimensional  projection  of  V. 

4Although  it  is  classically  referred  to  as  the  "Burgers"  equation,  the  first  time  that  it  was  introduced  dates  back 
to  the  1915  work  by  Bateman  [4], 
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(16) 


To  proceed,  let  us  recast  Eqs.  (3)  in  compact  form  as 


<9q 

dt 


+  £(q)  =  “g 


where  q  =  [p,Ui,0]T  is  the  array  of  the  solution  variables  [i  —  1,2,3)  and  £(q)  contains  all  the 
differential  operators  that  are  easily  identifiable  in  the  system  above.  From  now  on,  the  symbols 
7  and  ~  that  identify  the  large  scale  quantities  will  be  omitted  to  relax  the  notation.  Within  each 
element,  the  finite  dimensional  approximation  of  q(x,f)  is  given  by  the  expansion 


(iV+l)3 

qfe(x,t)|ne=  Y  V’fe(x)qfe(t)|ne 

k= 1 


(17) 


where  (N  + 1)3  is  the  number  of  collocation  points  within  the  element  of  order  N  and  il>k  are  the 
interpolation  polynomials  evaluated  at  point  k.  The  basis  functions  t/’fe  are  constructed  as  the  tensor 
product  of  the  one-dimensional  functions  ft.Q(£(x)),  hp(r](x)),  and  ft7(<C(x))  as: 

V’fc  =  M£(x))®Mf?(x))®MC(x)),  Va,/3,7  =  0, 

The  functions  ha(£(x)),  hp(r](x)),  /i7( C(x))  are  the  basis  functions  associated  with  the  IV +  1 
Legendre-Gauss-Lobatto  (LGL)  points  £Q,  pp,  £7,  respectively,  which  are  given  by  the  roots  of 


(1-C2)4(«  =  0, 

where  P'n(£)  are  the  derivatives  of  the  ./V  -order  Legendre  polynomial  [14;  32],  Given  these  defi¬ 
nitions,  the  one-dimensional  Lagrange  polynomials  ha  p  y(£)  are 


h'Ct,P,'y{£)  — 


a -awo 


N(N+1)  (£-£Q,/3,7)-Pjv(£a,/3,7) 


The  same  expansion  (17)  is  used  to  construct  the  derivatives.  By  differentiation  of  Eq.  (17) 
with  respect  to  time  and  space,  we  write: 


<9q  fe(x,f), 
dt  1 

9qh(x,t) 
dx  1 


( N+iy 


=  Y  ^(x) 


k= 1 
(N+ 1)3 


dqfc(t) 

dt 


=  £ 


k= 1 


dx 


(18) 


The  definitions  and  expansions  introduced  above  yield  the  weak,  element  problem 


dnhe 


(19) 


By  virtue  of  the  global  assembly  procedure,  the  global  problem  is  solved  on  Qh  and  consists  in 
finding  q h(t,x)  G  Vh  such  that 


>C(qh) 


dQn  = 


!nh 


iphgdtlh  \/iph  G  Vh. 


(20) 


Notice  that  we  have  not  explicitly  treated  the  second  order  operators  so  far.  Let  us  start  with  the 
momentum  equation  and  treat  the  strain 


l3rfS  n  d  (du 


p  dxj 


p  dxj 


dxj 


duj_ 

dx; 


where  the  dissipative  coefficient  p  is  taken  as  a  local  constant  (at  the  current  time  step)  that  can 
hence  be  taken  out  of  the  derivation.  Multiplication  by  the  test  function  xf>h  and  integration  yields 


jhp  d  f  dm  du 
^*--5—  brr  +  7rr  )dn 


p dxj  \  dx 


dxj 


that  becomes,  after  integration  by  parts, 


dm  diij 


d(iph/p)  /  dui  duj 


dx.j 


\  dxj  dxi 


dilp 


(21a) 


(21b) 


where  ny,re  is  the  j  component  of  the  unit  normal  vector  on  the  element  boundary  and 


d(iph/p)  ldiph  hd(l/p) 

dxj  pdxj+P  dxj  '  1  J 

Because  the  stresses  in  the  current  method  represent  an  artificial  diffusion  that  is  only  an  approx¬ 
imation  to  the  actual  Navier-Stokes  stresses,  to  limit  the  operation  count  in  our  implementation 
we  omit  the  second  term  on  the  right  hand  of  side  of  (22).  Furthermore,  because  of  the  continuity 
of  the  solution  across  the  element  boundaries  and  because  of  the  inviscid  or  periodic  boundary 
conditions  that  we  will  apply  in  the  tests  below,  the  flux  term  that  arises  from  the  integration  by 
parts  vanishes. 

The  final  semi-discrete  matrix  system  results  from  the  global  integral  equation  (20)  using  the 
expansions  (18)  and  the  global  assembly  of  the  element  matrices  that  result  from  it.  The  system 
written  in  compact  matrix- vector  form  is 

^  =  M-1£fe(qfe),  (23) 

where  M  is  the  global  mass  matrix  that  is  diagonal  by  construction  since  the  integration  and 
interpolation  points  within  the  elements  are  co-located.  For  more  details  on  the  construction  of 
this  system  and  of  the  global  assembly  operation  see,  e.g.,  [33]. 

System  (23)  is  discretized  in  time  by  the  implicit-explicit  (IMEX)  time  integration  described 
in  [19].  Leaving  the  details  of  IMEX  to  [19],  the  IMEX  time  approximation  uses  an  implicit 
approximation  of  the  (linear)  terms  responsible  for  the  fast  moving  acoustic  and  gravity  waves, 
whereas  the  slow  non-linear  advection  is  treated  explicitly.  In  this  paper,  IMEX  is  based  on  a 
fourth  order  multi-stage  additive  Runge-Kutta  (ARK4)  method.  This  approach  allows  us  to  relax 
the  constraint  that  fast  waves  put  on  the  Courant-Friedrichs-Lewy  [12]  condition  and  on  the  size 
of  the  time-step. 


4  Model  verification:  Euler  equations 

In  the  following  section,  Dyn-SGS  is  tested  against  a  suite  of  standard  benchmarks  of  ubiquitous 
use  when  testing  new  atmospheric  dynamical  cores. 
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First,  we  perturb  a  neutrally  stable  atmosphere  with  thermal  anomalies  that  vary  in  definition 
and  size.  The  tests  are  divided  into  2gD,  where  the  domain  extends  to  infinity  along  the  y-direction, 
and  fully  3D  problems  in  a  simply  connected  domain.  They  include  the  following:  (i)  2^D  rising 
thermal  bubble  in  a  large  domain  (see  the  original  2D  test  in  [2])  and  (ii)  the  classical  density 
current  [48].  In  3D  we  have:  (iii)  a  rising  thermal  bubble  in  a  small  domain  [33],  (iv)  the  baroclinic 
wave  in  a  channel  [51].  Except  for  the  balanced  initial  state  of  the  baroclinic  wave  in  a  channel, 
an  analytic  solution  does  not  exist  for  these  problems.  For  this  reason,  it  must  be  understood  that 
most  of  these  tests  can  only  give  a  qualitative  (and  relative)  information  on  the  accuracy  that  one 
model  can  achieve. 

4.1  2^D  rising  thermal  bubble  in  a  large  domain 

This  test  consists  of  a  flow  that  is  triggered  by  the  thermal  perturbation  of  a  neutrally  stratified 
atmosphere  at  initially  uniform  potential  temperature  6q  =  300  K  and  in  hydrostatic  equilibrium 
such  that  the  pressure  decreases  with  2  as: 


P  =  P0 


Cp  /  R 


(24) 


The  domain  Q  =  [—5000,5000]  x  [—00,00]  x  [0, 10000]  m3  and  the  definition  of  the  perturbation  are 
given  as  in  [2] .  The  perturbation  is  linear  and  defined  as 


a  e  =  ec 


if  r  <  ro  =  2000m, 


(25) 


where  r  —yj (x  —  xc )2  +  (z  —  zc)2,  (xc,zc)  =  (0,2000)m,  and  9C  =  2  K.  Due  to  the  symmetry  of  the 
problem,  we  simulate  only  half  of  the  domain  (to  verify  that  the  current  method  can  indeed  preserve 
symmetry,  in  Section  4.3  we  analyze  a  fully  3D  simulation  without  any  geometric  assumptions). 
The  initial  velocity  field  is  zero  everywhere.  Periodic  boundary  conditions  are  used  along  y  whereas 
no-flux  is  imposed  in  x  and  z.  We  perform  four  runs  on  four  different  grids  with  effective  resolutions 
1)  Ax  =  Az  =  125m  (as  in  [2]),  2)  Ax  =  A z  =  62.5m,  3)  Ax  =  Az  =  31.25m,  and  4)  Ax  =  Az  = 
15.625m.  The  different  resolutions  were  used  to  analyze  the  behavior  of  the  method  as  the  grid 
is  refined,  although  no  proper  convergence  study  is  made.  The  contour  lines  of  the  perturbation 
potential  temperature  A 9  are  plotted  at  t  =  1020  s  in  Fig.  1.  The  maximum  of  A 9  at  the  final 
time  is  ss  1.4  K,  which  agrees  with  the  /-wave  solution  of  [2]  at  the  same  resolution  of  125  m,  and 
is  equivalent  with  the  WRF  simulation  reported  in  the  same  paper.  The  extrema  of  A 9  in  this 
study  increase  a  few  fractions  of  a  degree  as  the  grid  is  refined  to  15.625  m.  This  is  to  be  attributed 
to  the  grid  dependence  of  the  current  SGS  model,  whose  dissipation  properties  vary  quadratically 
with  respect  to  the  size  of  the  grid.  We  will  touch  more  on  this  issue  in  the  analysis  of  the  density 
current  below.  To  maintain  the  WRF  solution  stable,  Ahmad  and  Lindeman  had  to  add  a  constant 
15m2  s_1  diffusion.  They  did  not  for  their  finite  volume-based  /-wave  decomposition.  Without 
a  constant  diffusion,  we  preserve  stability  for  all  the  tested  resolutions.  The  shear  between  the 
boundary  of  the  perturbation  and  the  background  at  rest  triggers  very  well  defined  Rayleigh- Taylor 
instabilities  that  become  ever  more  visible  when  the  resolution  is  increased.  To  give  some  hints  to 
the  reader  on  what  the  sub  grid  scales  look  like  and  what  spatial  distribution  they  have  with  respect 
to  the  developing  bubble,  we  show  them  in  Fig.  2.  In  the  figure,  the  symbol  SGS(-)  indicates  the 
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stresses  that  appear  in  the  momentum  and  9  equations.  They  are  (with  some  abuse  of  notation) 

drSGS  80 

SGS(ui)  =  -jj—/IISGS(ui)ll00in,  and  SGS(9)  =  -^/||SGS(0)||oo>n. 

OXj  OXj 

In  [43]  the  Prandtl  number  is  assigned  the  artificial  value  0.1.  Because  this  value  is  fully  artificial 
and  does  not  have  a  physical  foundation,  for  this  study  we  decided  to  test  how  the  current  stresses 
behave  if  we  picked  the  air  Prandtl  number,  approximately  equal  to  0.7.  We  also  tried  0.1  and  an 
intermediate  value  of  0.35.  With  Pr=0.1  the  solution  still  preserves  stability  for  all  the  flow  regimes 
that  we  tested,  confirming  that  the  method  was  implemented  properly  since  0.1  was  successfully 
used  in  [43].  However,  in  our  opinion  the  solution  lost  some  of  its  physics  since  hardly  no  dissipation 
of  A 9  could  be  observed  as  time  evolved.  In  other  words,  the  maximum  of  A 9  was  preserved  along 
the  full  simulation,  which  is  an  unphysical  behavior  for  the  flows  we  are  interested  in  simulating. 
If  we  looked  at  this  from  a  mathematical  point  of  view,  this  is  arguably  a  great  result  since  the 
method  achieved  stability  without  really  affecting  the  originally  inviscicl  system;  however,  from  a 
more  physical  point  of  view,  we  want  a  method  that  not  only  stabilizes,  but  also  introduces  the 
necessary  sub-grid  stresses  that  model  the  physical  dissipation  of  a  moving  atmosphere  in  a  sheared 
environment.  In  Fig.  3  we  compare  the  solution  of  the  same  rising  thermal  bubble  already  shown 
in  Fig.  1  for  Pr=0.35,  but  now  also  using  Pr=0.7  (top  row  of  Fig.  3).  As  expected,  Pr=0.7  gave 
a  smoother  solution,  possibly  indicating  a  greater  dissipation.  However,  excluding  the  extremely 
small  oscillation  that  are  visible  in  certain  regions  in  the  case  of  Pr=0.35,  we  observe  how  the 
distributions  of  A 9  for  the  two  cases  are  practically  identical.  The  extrema  of  A 9  are  preserved; 
this  is  seen  by  looking  at  the  differences:  [A0maXipr=o. 35  —  A9maX}pr=0. 7]  (15.625m)  =  0.028  K  and 
[A0maa,!pI.=o.35  —  A0maX!pr=o. 7] (31.25m)  =0.1  K;  they  possibly  show  that  the  more  physical  Pr 
and  the  dynamically  adapting  diffusion  coefficients  have  an  effect  on  the  smaller  scales  (i.e.  small 
oscillation  in  the  domain)  rather  than  on  the  global  solution  and  indeed  help  smoothen  the  solution 
only  where  oscillations  are  important. 

4.2  2^D  Density  current 

The  density  current  was  introduced  in  [10]  and  became  a  standard  benchmark  in  the  development 
of  atmospheric  codes  [48] .  Like  in  [2] ,  in  this  paper  the  benchmark  is  run  without  the  constant  and 
uniform  artificial  diffusion  with  coefficient  /.t  =  75m2s~1  of  [48].  This  is  because  we  are  interested 
in  assessing  Dyn-SGS  as  a  stabilizing  tool  that  does  not  require  additional  constant  viscosity.  For 
this  reason,  no  converged  solution  should  be  expected.  On  the  contrary,  as  the  grid  is  refined,  more 
structures  will  be  resolved.  We  will  show  this  shortly.  The  background  initial  state  is  at  a  uniform 
potential  temperature  9q  =  300  K  within  the  domain  H  =  [—25600,25600]  x  [—00,00]  x  [0,6400]m3. 
A  perturbation  of  9  centered  in  (xc,zc)  =  (0,3000)m  and  with  radii  (rx,rz)  =  (4000, 2000) m  is  given 
by  the  function 

Q 

A9  =  -A  [1  +  cos(7rcr)]  if  r  <  1  (26) 

where  9C  —  — 15  K  and  r  =\J [x  —  xc)/rx  +  (z  —  zc) / r 2.  Periodic  boundary  conditions  are  used  along 
y  whereas  no-flux  and  free-slip  conditions  are  set  in  x  and  z.  The  initial  velocity  is  zero.  The 
time  evolution  of  the  density  current  for  5400  s  is  shown  in  Fig.  4.  Because  A 9  is  transported 
(and  diffused)  by  the  flow,  it  is  the  perfect  variable  to  be  used  in  the  visualization  of  mixing.  The 
evolution  and  transition  onto  a  fully  mixed  flow  is  evident.  The  mixing  is  triggered  by  the  shear 


11 


Z  (m)  Z  (m) 


12 


10000 


7500 


5000 


2500 


10000 


7500 


5000 


2500- 


Current 
t=1020  s 


31.25  m 

2500 

X  (m) 


N 


10000 


7500 


5000 


2500- 


N 


10000 


7500 


5000 


2500 


5000 


5000 


Current 
t=1020  s 


15.625  m 

2500 

X  (m) 


5000 


Figure  1:  2^D  rising  thermal  bubble  in  a  large  domain.  Four  resolutions  as  indicated  in  each 
subfigure.  The  perturbation  of  potential  temperature  Ad  is  plotted  at  t  =  1020  s.  All  the  results 
are  shown  for  a  maximum  Courant  number  8. 
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Figure  2:  2|D  rising  thermal  bubble  in  a  large  domain.  From  left  to  right:  contours  of  the  dynamic 
stabilization  stresses  SGS(u ),  SGS(w ),  and  SGS(d).  The  values  of  the  SGS  terms  are  normalized 
with  respect  to  the  global  maximum  of  the  stress  at  hand  (i.e.  SGS (9)  is  the  stress  of  the  9  equation 
and  is  normalized  with  respect  to  its  maximum  value  at  the  given  time  step).  These  plots  show 
the  structure  of  these  stresses  and  how  they  dynamically  adapt  to  the  solution.  These  stresses  are 
obtained  from  the  fine  grid  simulation  using  the  resolution  Ax  =  Az  =  15.625  m. 


that  becomes  greater  as  the  velocity  of  the  moving  cold  mass  of  air  first  accelerates  and  later  hits 
the  right  wall.  It  bounces  off  the  boundary  and  keeps  moving  left-  and  rightwards  as  more  mixing 
occurs  at  the  expense  of  kinetic  energy.  Due  to  mixing  and  dissipation,  the  minimum  of  A 9  increases 
from  -15  K  at  t  =  0  s  to  -3.729  K  during  5400  seconds  (from  dark  blue  to  light  green  shading  in  the 
plots),  whereas  the  maximum  is  preserved  to  0  K  (dark  red). 

Figure  5  shows  A 9  for  the  effective  resolutions  Ax  =  Az  =  50  m,  Ax  =  Az  =  25  m,  and  Ax  = 
Az  =  12.5  m  at  simulation  times  600,  750,  and  900  s.  To  our  knowledge,  no  data  is  available  to 
compare  against  at  the  12.5  m  resolution.  It  is  evident  that  the  amount  of  vortical  structures  is 
larger  at  higher  resolution.  Certainly,  the  finer  resolution  is  expected  to  yield  a  more  resolved 
solution.  However,  we  believe  this  to  be  a  mixed  effect  of  the  resolution  and  smaller  dissipation 
required.  Just  like  the  classical  Smagorinsky  and  other  models,  the  current  SGS  model  is  a  quadratic 
function  of  the  grid  size  so  that  dissipation  decreases  as  the  grid  is  refined.  The  effect  on  the  solution 
is  non-linear.  The  non-linearity  can  be  explained  as  follows:  the  finer  grid  lowers  the  influence  of 
dissipation  with  a  direct  effect  on  the  oscillations  of  the  solution  variable.  Slightly  larger,  yet 
controlled,  localized  oscillations  imply  larger  gradients  (and  hence  residuals)  that  directly  affect  the 
dissipation  in  the  regions  of  larger  residuals. 

Given  the  accepted  use  of  the  Smagorinsky  model  [37;  47]  in  numerical  weather  prediction,  to 
support  our  hypothesis  that  a  properly  designed  SGS  model  can  serve  as  a  stabilization  method, 
we  compare  the  results  of  the  current  model  with  the  constant  coefficient  Smagorinsky  at  25  m  and 
50  m  resolutions.  The  results  are  plotted  in  Figs.  6-7.  The  two  methods  are  in  strong  agreement 
at  both  resolutions.  It  must  be  kept  in  mind  that  the  density  current  described  here  is  not  a 
boundary  layer  flow;  the  boundaries  are  treated  as  if  the  problem  were  inviscid  although  the  flow 
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Figure  3:  As  Fig.  1,  but  here  we  concentrate  on  the  effect  of  the  value  of  the  Prandtl  number.  Top 
row:  Pr=0.7.  Bottom  row:  Pr=0.37.  Left  column:  Ax  =  Az  =  31.25  m.  Right  column:  Ax  =  Az  = 
15.625  m.  As  expected  the  Pr=0.7  solution  is  smoother,  possibly  indicating  a  greater  dissipation. 
However  the  differences  [A9rnax^pr= 0.35  ^@max,Pr= 0.7] ( 15.625m)  0.028  K  and  [A9rnax,pr= 0.35 

A0maXjpr=o. 7]  (31.25m)  =  0.1  K  show  that  the  greater  Pr  has  an  effect  on  the  smaller  scales  only 
rather  than  on  the  global  solution. 
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Figure  6:  2^D  density  current.  Comparison  of  the  current  method  (top)  with  the  solution  obtained 
using  Smagorinsky  with  constant  Cs  =  0.14  (bottom)  using  a  grid  resolution  Aa:  =  A z  =  50  m. 


is  viscous;  free-slip  boundary  conditions  are  applied  on  every  solid  wall.  This  is  an  important  point 
since  it  is  well  known  that  Smagorinsky  is  overly  dissipative  in  boundary  layer  flows  unless  it  is 
properly  corrected.  We  are  currently  working  on  a  thorough  comparison  of  three  LES  models  that 
are  commonly  used  in  atmospheric  simulations  and  the  results  will  be  reported  in  a  subsequent 
paper. 

This  test  case  does  not  admit  an  analytic  solution.  To  compare  with  other  models  we  look 
at  the  front  position  at  t  =  900  s,  where  the  front  is  defined  as  the  position  on  the  ground  where 
A 8  =  —1  K.  In  Table  3,  the  position  of  the  front  is  reported  for  Dyn-SGS,  for  our  implementation 
of  Smagorinsky,  and  for  some  of  the  results  reported  in  the  literature.  As  the  grid  is  coarsened,  the 
front  is  slightly  slower;  this  fact  is  also  observed  in  Fig.  5  of  [48].  On  the  50  m  grid,  all  models 
agree  with  a  front  location  in  the  range  [14409, 14975]  m.  For  the  data  that  are  available  at  different 
resolutions  (see  current,  Smagorinsky,  and  VMS  in  the  table),  a  trend  is  observed:  as  the  resolution 
is  increased  from  200  m  to  12.5  m,  the  front  moves  relatively  faster  although  within  a  few  meters 
of  difference  in  the  front  position  from  the  fine  to  the  coarse  grid  solution.  The  somewhat  smaller 
speed  in  the  case  of  a  coarser  grid  is  to  be  attributed  to  the  relatively  larger  dissipation  in  the  coarse 
solution;  this  makes  the  flow  more  viscous  and  hence  gives  it  a  tendency  to  be  slightly  slower. 

As  the  resolution  is  increased,  the  amount  of  structures  that  are  resolved  increases  as  well. 
Without  the  large  viscosity  that  homogenizes  the  solution  as  done  in  [48]  with  the  sole  target  of 
reaching  convergence,  the  inviscid,  non-linear,  and  non-steady  solution  that  we  present  here  is  not 
expected  to  show  signs  of  space-convergence.  The  same  behavior  was  observed  in  [41],  where  VMS 
was  used  to  stabilize  a  finite  element  discretization  of  the  Euler  equations.  Rather,  we  expect  more 
and  more  structures  to  be  resolved  until  a  grid  resolution  of  the  order  of  the  smallest  eddies  is 
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Figure  7:  As  in  Fig.  6,  but  at  the  finer  resolution  Ax  =  A z  =  25  m. 


reached. 


4.3  3D  warm  bubble 

So  far  only  2^D  problems  have  been  presented,  where  the  solution  was  sought  in  a  3D  domain  with 
infinite  extension  along  the  y-direction.  To  test  the  complete  applicability  of  the  current  method 
to  fully  three-dimensional  problems,  we  consider  the  3D  analog  of  the  2^D  buoyant  rising  thermal 
described  above.  The  problem  is  now  solved  in  the  domain  D  =  [0,1000pm3  as  in  [33].  The  initial 
perturbation  Ad  is  no  longer  linear;  it  is  given  by  the  function 


1.0  + cos 


if  r  <  ro  =  250  m, 


(27) 


where  r  =y/ ( x  —  xc)2  +  {y  —  yc )2  +  (z  —  zc)2,  ( xc,yc,zc )  =  (500,500,260)  m,  and  9C  =  0.5  K.  To  verify 
the  preservation  of  the  axial  symmetry  of  the  problem,  we  solve  the  fully  3D  problem  without  relying 
on  its  axial  symmetry.  In  Figs.  8  and  9,  the  4th  and,  respectively,  8  -order  solutions  are  shown 
at  t  =  400  s.  With  ne  =  10  x  10  x  10,  the  effective  resolutions  are  25  and  12.5  m.  Based  on  [33], 
a  constant  and  uniform  artificial  diffusion  (indicated  by  the  symbol  HV 2)  is  expected  to  lead  to 
a  solution  that  varies  according  to  the  value  of  its  coefficient.  The  current  method  is  designed  to 
dissipate  the  solution  only  where  necessary.  The  value  of  the  extrema  of  the  current  solution  is 
in  close  agreement  with  the  HV 2  solution  using  0.5m2s_1  (max(A0)  «  0.35),  indicating  that  the 
amount  of  diffusion  that  is  being  added  for  stabilization  is  sufficiently  small  to  not  dissipate  the 
solution  away  from  the  region  of  larger  residuals,  but  is  large  enough  to  suppress  all  the  oscillations 
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Table  3:  Density  current.  Comparative  results  of  front  location  at  t  =  900  s.  The  results  are 
reported  for  the  following  models:  current  model,  Smagorinsky  implemented  in  our  code,  VMS, 
WRF-ARW,  /-wave,  filtered  Spectral  Elements  (SE),  filtered  Discontinuous  Galerkin  (DG),  REFC, 
REFQ,  and  PPM.  All  models  but  the  current,  Lilly-Smagorinsky,  and  VMS,  used  artificial  diffusion 
with  constant  12  =  75  m2  s_  1 . 


Model 

Space  discr. 

Resolution 

ne 

Order 

/i  =  75nVs  1 

Front  Location  [m] 

Current  (Dyn-SGS) 

SEM 

12.5  m 

512  x  lx  128 

4  th 

No 

15056 

" 

" 

25  m 

256  x  1  x  64 

^th 

No 

14992 

11 

" 

50  m 

128  x  1  x  32 

^th 

No 

14535 

" 

" 

100  m 

64  x  1  x  16 

^th 

No 

14325 

" 

" 

200  m 

32  x  1  x  8 

^th 

No 

13552 

11 

" 

133  m 

32  x  1  x  8 

6th 

No 

14568 

" 

" 

100  m 

32  x  1  x  8 

8th 

No 

14754 

Smagorinsky 

SEM 

25  m 

256  x  1  x  64 

4  th 

No 

14918 

11 

" 

50  m 

128  x  1  x  32 

^th 

No 

14726 

" 

" 

100  m 

64  x  1  x  16 

^th 

No 

14551 

VMS  [41] 

FEM 

25  m 

No 

14890 

"  [41] 

" 

50  m 

No 

14629 

"  [41] 

" 

75  m 

No 

14487 

"  [41] 

" 

100  m 

No 

14355 

SE  [20] 

SEM 

50  m 

Yes 

14767 

DG  [20] 

DG 

50  m 

Yes 

14767 

/-wave  [2] 

FV 

50  m 

Yes 

14975 

WRF-ARW  [2,  run  in-] 

FD 

50  m 

Yes 

14470 

REFC  [48] 

FD 

50  111 

Yes 

14437 

REFQ  [48] 

FD 

50  m 

Yes 

14409 

PPM  [48] 

FD 

50  m 

Yes 

15027 

that  are  visible  in  the  case  of  HVi-  As  previously  observed  for  the  density  current,  the  higher  the 
resolution  the  faster  the  front.  Although  the  difference  in  height  is  minimal,  it  is  still  noticeable 
within  a  few  meters  difference.  This  same  test  was  run  in  [33]  using  an  effective  resolution  of  12.5  m 
(8th-order)  and  in  [19]  with  an  effective  resolution  of  10  m.  In  the  latter,  the  front  is  approximately 
50  m  higher  than  it  is  in  the  case  of  [19].  This  said,  the  current  method  produces  a  thermal  whose 
dissipation  is  as  limited  and  constrained  in  space  as  it  is  for  I/V2(0.5m2s_1)  but,  at  the  same 
time,  preserves  the  smoothness  of  FfV2(5.0m2s^1).  The  4th  and  8  -order  solutions  show  the  same 
behavior.  It  is  also  shown  that  its  stability  properties  are  not  compromised  or  modified  as  the  order 
of  the  spectral  element  is  increased  from  4  to  8.  It  is  important  to  keep  this  in  mind  since  the 
increase  of  the  element  order  coincides  with  an  important  decrease  in  the  grid  spacing  within  each 
element,  and  the  current  version  of  the  stabilizing  method  is  a  function  of  the  element  effective 
resolution. 

4.4  Baroclinic-wave  in  a  3D  channel  with  geostrophically  balanced  back¬ 
ground 

A  geostrophically  balanced  background  is  defined  as  in  [51].  The  /—plane  approximation  is  consid¬ 
ered.  The  flow  is  confined  in  a  very-high  aspect  ratio  domain  fl  =  [Lx  x  Ly  x  Lz\  =  [40000  x  6000  x  30]  km3 
with  periodic  boundary  conditions  along  the  flow  direction  and  free-slip,  non-penetrating  bound¬ 
aries  everywhere  else.  Notice  how  the  characteristic  size  of  this  domain  (and  of  the  grid  spacing) 
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Figure  8:  3D  rising  thermal  bubble.  Solution  at  t  =  400  s.  Left  column:  irz-slice  of  Ad  at  y  =  500 
m.  Right  column:  yz-slice  at  x  =  500  m.  Top  row:  solution  using  the  current  method  (Dyn-SGS). 
Middle  row:  solution  using  second-order  hyper  diffusion  (HV2)  with  constant  coefficient  /.<  =  0.5. 
Bottom  row:  second-order  hyper  diffusion  with  constant  coefficient  (i  =  5.0.  Courant  number  20.0 
on  a  25  m  effective  resolution  grid  of  order  4.  Q  =  [1000] 3  m3. 
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Figure  9:  As  Fig.  8,  but  for  8  -order  elements  and  effective  resolution  of  12.5  m,  and  Courant 
number  17.5. 


is  up  to  0(1O4)  times  larger  than  that  of  the  previous  problems.  If  we  consider  a  pressure-based 
vertical  coordinate  77,  from  which  the  z  coordinate  is  derived  via  iteration  (see  appendix  in  [30]), 
the  initial  jet  is  a  zonally  symmetric  flow  defined  as 

u(x,y,r j)  =  —  itosin2  (  — ■ —  )  ln^exp 

V  Ly  J 


with  amplitude  uq  =  35 ms^1  and  vertical  width  parameter  b  =  2.  The  meridional  and  vertical 
velocities,  (v,w),  are  initially  zero.  In  the  /—plane  approximation  the  Coriolis  parameter  is  fo  = 
2wsin</?o  at  the  latitude  ipo  =  45°N,  where  oj  is  the  Earth  rotational  velocity  (Table  1).  Although 
in  a  plane  channel,  this  jet  is  designed  to  resemble  a  mid-latitude  westerly  zonal  wind  with  a  zonal- 
and  time-mean  jet  speed  at  the  earth  troposphere.  The  background  geopotential  is  given  by 


^(x,y,rf) 


^1  —  yRdr/9^j  +  A$(x,y)ln77exp 


(29) 


where  Tq  =  288  K  is  a  reference  temperature,  T  =  0.005 Km  1  is  the  lapse  rate,  and  Ad*  is  the 
perturbation  of  $  given  by 


Ad>  = 


The  perturbed  temperature  distribution  is  given  by 


T{x,y,rj)  =  T0r]Rdr/9  +  A<^’  ^  exp 


(30) 


(31) 


Baroclinic  instabilities  are  responsible  for  mid-latitude  cyclones  [26]  and  are  thus  important  atmo¬ 
spheric  processes  for  an  atmospheric  model  to  capture.  The  baroclinic  wave  instability  is  triggered 
by  a  perturbation  of  the  initially  balanced  zonal  velocity  field.  As  the  wave  breaks,  gravity  waves 
are  radiated  with  the  intent  of  restoring  the  initial  geostrophic  balance  [25].  The  perturbation  is 
given  by  an  unbalanced  smooth  profile  centered  at  (xc,yc)  =  (2000, 2500)  km  and  defined  as 


A u(x,y,z)  =  Up  ex p 


(x  ~  xc)2  +  (y  -  yc)2 
Lp 


(32) 


where  up  =  lms_1  is  the  perturbation  amplitude  and  Lp  =  600  km  is  the  width  parameter. 

We  have  seen  above  how  the  dynamic  dissipation  depends  on  a  characteristic  grid  size.  Due  to 
the  very  high  aspect  ratio  of  the  current  grids,  with  the  vertical  resolution  being  100  times  smaller 
than  the  horizontal  one,  stabilization  for  this  problem  was  first  run  by  setting  the  vertical  diffusion 
to  a  value  100  times  smaller  along  2  than  it  is  along  x  and  y ,  and  then  simply  to  zero  along  z.  The 
result  did  not  change  significantly  so  that  the  plots  shown  in  this  paper  are  only  those  obtained  with 
zero  vertical  dissipation.  The  vertical  resolution  is  sufficiently  high  to  preserve  stability  without 
the  need  for  vertical  diffusion. 

To  get  a  sense  of  the  error  that  we  commit  and  compare  against  [51],  we  first  ran  this  test 
without  perturbing  the  initially  balanced  flow.  The  relative  error  norms  Li,  L2,  and  of  q =p9 

are  defined  as 


Li  (<z(x,t)) 


/nfa|g(x,f)-g(0,x)|dflfe 

fQh  \q(0,x)\dnh 
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Table  4:  Error  norms  of  p9  in  the  solution  of  the  geostrophically  balanced  flow  in  a  3D  channel  for 
5  different  horizontal  resolutions,  but  the  same  1  km  vertical  grid  spacing.  The  errors  are  computed 
at  day  1  with  respect  to  the  initial  condition.  For  comparison,  the  errors  of  the  non- viscous  solution 
are  also  reported  and  are  indicated  by  the  ^  symbol. 


km 

Lx 

Lx 

l2 

l2 

Loo 

Loo 

800 

2.77x  10_® 

2.75x10“® 

7.44xl0_lu 

7.26x  10_iu 

4.27x10"® 

4.19x10“® 

400 

3.11x10-® 

3.11x10-® 

1.14X10”11 

1.14x  10-11 

6.91x10-® 

6.91x10-® 

200 

8.81x10-® 

8.81x10-® 

6.27xl0~15 

9.27x  10-15 

4.03x  10-7 

4.03x  10-7 

100 

7.21X10"9 

7.21xl0“9 

5.35xl0~17 

5.35xl0~17 

1.78x10-® 

1.78x10-® 

50 

6.26xlO-10 

6.25xlO-10 

4.16X10-19 

4.16x  10-19 

1.24x  10-9 

1.24x  10-9 

and 


L2  (g(x,f))  = 


1  fnh(q(x,t)-q(0,x))  dflh 
fQhq2(0,x)dnh 


Loo  (q(x,t)) 


max \q(x,t)  —  q{  0,x)  | 
max  |^(0,  x)  | 


where  q( 0,x)  is  the  initial  solution  and  the  integrals  are  computed  via  the  usual  quadrature  formulas. 
We  report  their  values  at  day  1  in  Table  4.  The  errors  obtained  with  stabilization  are  compared 
against  our  inviscicl  solution  and  against  the  inviscid  solution  of  [51].  Because  this  problem  is 
smooth  and  the  flow  is  stationary,  the  dissipative  effect  of  the  current  method  should  be  negligible. 
Although  some  small  dissipation  is  still  partially  active  across  the  whole  domain  (plot  of  the  SGS 
is  not  shown  for  this  case),  its  effect  is  indeed  minor,  as  can  be  seen  by  the  almost  non-existent 
difference  between  the  tabulated  errors  of  the  stabilized  and  inviscid  solutions.  These  errors  are 
plotted  in  Fig.  10  as  well;  the  error  decay  shows  that  the  order  of  the  numerical  approximation  is 
not  degraded.  In  the  same  figure,  the  time  evolution  of  the  error  is  plotted  for  a  20  day  simulation. 

In  the  case  of  the  perturbed  flow  that  triggers  a  baroclinic  wave,  the  evolution  of  T  from  day  12 
to  day  14  is  shown  in  Fig.  11  given  an  anisotropic  grid  with  Ax  =  100  km,  Ay  =  75  km,  Ax:  =  1.25 
km.  (Note:  to  our  knowledge,  no  result  has  been  shown  past  day  12.  We  believe  an  instability  in  the 
problem  statement  is  responsible.)  The  plots  show  the  solution  on  the  xy  cross  section  at  z  =  500 
m.  With  respect  to  the  results  of  [51],  the  wave  breaking  occurs  at  the  same  time  (approximately 
9  days).  The  values  of  T  and  of  the  vertical  component  of  vorticity  are  in  agreement  with  [51]  at 
all  resolutions.  The  flow  at  t  >  12  days  is  shown  as  well  for  future  comparisons. 


4.5  Mass  loss 


Next,  we  wish  to  analyze  the  effect  of  Dyn-SGS  on  the  mass  conservation  since  this  is  an  important 
metric  for  atmospheric  models.  Let  us  define  the  time  dependent  normalized  mass  loss  as 


M{t)ioss 


JnP(to)dtt 


(33) 


where  pit o)  indicates  density  at  the  initial  time.  Figure  12  shows  M(t)ioss  for  a  20400  s  simulation 
of  the  2^D  rising  thermal  bubble  at  three  different  resolutions  and  three  stabilization  methods:  the 
current  one,  a  2nd-orcler  artificial  diffusion,  HV2  with  constant  coefficient,  and  a  4t,l-order  hyper 
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Figure  10:  Normalized  h-e rror  norms  for  the  solution  of  the  geostrophically  balanced  flow  in  a  3D 
channel.  Left  plot:  grid  convergence  (log-log  scale).  Right:  time  evolution  of  1 1 1 1 2 - 


diffusion,  HV4,  with  constant  coefficient.  All  methods  show  a  trend  towards  higher  mass  loss  values 
as  the  resolution  is  increased.  For  instance,  the  current  case  has  the  lowest  value  that  approximates 
machine  precision  in  the  case  of  85  and  170  m  resolutions  but  increases  to  higher  values  in  the  case 
of  the  45  m  resolution.  In  the  case  of  HV 2  and  HV4,  the  mass  loss  is  always  in  the  proximity  of 
1  x  1CT14,  although  the  curves  are  visibly  higher  for  higher  resolutions.  In  spite  of  the  differences 
among  the  three  methods,  a  sustained  mass  loss  0(<  10~14)  is  an  acceptable  value.  At  the  coarsest 
resolution  of  170  m,  HV4  failed  to  preserve  stability  with  the  given  coefficient;  this  is  why  the 
time  series  is  truncated  at  approximately  16000  seconds.  Although  a  different  coefficient  would 
immediately  solve  this  problem,  the  search  for  a  better  coefficient  falls  beyond  the  scope  of  this 
paper. 

In  Fig.  13,  we  plot  the  mass  loss  for  a  20  day  simulation  of  the  geostrophically  balanced  flow 
in  a  3D  channel.  The  very  large  time  scale  of  this  problem  is  a  good  test  for  mass  conservation 
in  a  geophysical  flow.  Although  not  at  machine  accuracy  and  in  spite  of  a  light  tendency  towards 
mass  increase  that  begins  at  approximately  day  6,  0(1  x  1CP15)  mass  loss  is  still  indicating  that 
the  method  is,  in  this  respect,  robust  for  very  long  simulations.  A  significant  mass  loss  increase  or 
oscillation  would  be  indicative  of  an  improper  behavior  of  the  model. 


5  Model  verification:  quasi-linear  scalar  equations 

We  test  the  model  against  the  solution  of  the  quasi-linear  ID  Burgers  equation.  The  reason  for 
including  this  test  is  twofold:  first,  we  demonstrate  that  Dyn-SGS  can  also  be  used  to  stabilize 
scalar  equations,  which  are  typically  required  as  part  of  atmospheric  modeling  systems;  second,  this 
problem  admits  an  exact  solution  to  compare  against.  The  solution  of  linear  transport  problems 
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Figure  11:  Baroclinic  instability  in  a  channel,  /-plane  solution  at  z  =  500  m.  The  temperature  is 
plotted  at  days  12,  13,  14  on  a  grid  resolution  Aa:  =  100  km,  Ay  =  75  km,  A z  =  1.25  km  using 
4 "''-order  polynomials.  The  thick  contour  corresponds  to  T  =  290  K. 
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Figure  12:  Time  evolution  of  the  mass  loss  for  the  stabilized  solution  using  the  current  method, 
using  HV 2  with  constant  coefficient  v  =  25m2 s_1,  and  using  HV \  with  v  =  l£+5m4s_1  in  the 
simulation  of  the  2^D  rising  thermal  bubble.  All  methods  show  a  trend  towards  higher  mass  loss 
values  as  the  resolution  is  increased.  For  instance,  the  current  case  has  the  lowest  value  that 
approximates  machine  precision  in  the  case  of  85  and  170  meter  resolutions,  but  increases  to  higher 
values  at  45  m  resolution.  In  the  case  of  HV 2  and  HV 4,  the  mass  loss  is  always  in  the  proximity 
of  1  x  10“ 14 ,  although  the  curves  are  visibly  higher  at  higher  resolution.  At  the  coarsest  resolution 
of  170  m,  HV4  failed  to  preserve  stability  with  the  given  coefficient;  this  is  why  the  time  series 
is  truncated  at  approximately  16000  seconds.  A  different  coefficient  would  immediately  solve  this 
problem,  although  its  search  is  not  relevant  in  the  context  of  this  paper. 
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Figure  13:  Time  evolution  of  the  mass  loss  for  the  stabilized  solution  of  the  geostrophically  balanced 
flow  in  a  3D  channel. 


via  Dyn-SGS  is  described  in  a  recent  report  [39]. 

The  solution  is  computed  in  il('x)  =  [0,2]  for  the  initial  condition  u(x,  0)  =  0.5 +  sin(7rx)  and 
periodic  boundary  conditions.  The  solution  is  advanced  in  time  until  t  =  2/7T.  The  solution  is 
computed  on  a  40  and  on  an  80  element  grid  of  order  8,  which  correspond  to  the  effective  resolutions 
Aa:  =  6.25  x  10~3  m  and  Aa’  =  3.125  x  10-3  m,  respectively.  The  exact  solution  of  the  inviscid 
Burgers  equation  is  computed  via  the  method  of  characteristics.  The  characteristic  curves  of  the 
computed  solution  are  plotted  in  Fig.  14.  As  time  approaches  l/ir  s,  the  characteristic  curves  begin 
to  cross,  generating  a  shock  that  becomes  stronger  as  time  passes.  At  all  times,  the  solution  remains 
smooth  sufficiently  far  from  the  shock,  as  visible  from  the  point-wise  relative  error,  \u  —  uexact\i 
plotted  in  Fig.  15  on  the  x  —  t  plane.  The  superposition  of  the  computed  and  exact  solutions  at 
the  two  time  levels  t  =  [3/(27t),  2/7t]  is  plotted  in  Figs.  16  and,  respectively,  17.  The  shock  is  well 
captured  within  one  element.  In  spite  of  the  high-order  approximation  (8  -order  elements),  the 
over-  and  under-shoots  in  the  proximity  of  the  shock  are  limited,  very  well  controlled,  and  do  not 
propagate  from  the  shock  across  the  domain.  This  is  confirmed  by  the  point-wise  error  curves  plotted 
in  Fig.  18.  Our  results  agree  with  the  entropy- viscosity  solution  of  [22].  Moreover,  the  sharpness  of 
the  discontinuous  solution  seems  uncompromised  by  the  action  of  the  dynamic  diffusion.  This  is  a 
strong  result  that  demonstrates  well  the  shock-capturing  capabilities  of  Dyn-SGS  when  high-order 
methods  are  used. 

Figure  19  shows  the  element-wise  structure  of  fi,  Eq.  (15),  at  t  =  [l/7r,  3/(27t),  2/7t]  s.  h  has 
no  effect  on  the  solution  in  the  smooth  regions,  as  it  clearly  acts  only  in  the  neighborhood  of  the 
shock. 

The  normalized  Li,L2,  and  Loo  error  norms  are  plotted  in  Fig.  20  as  a  function  of  time.  As 
expected,  the  formation  of  the  shock  has  a  major  effect  on  the  accuracy  of  the  solution,  as  is  visible 
from  the  jump  in  the  infinity  norm. 


6  Conclusions 

We  presented  the  application  of  a  dynamic  sub-grid  scale  model  (Dyn-SGS)  for  Large  Eddy  Sim¬ 
ulation  to  stabilize  the  spectral  element  (SEM)  solution  of  low  Mach  number  stratified  flows  and 
of  quasi- linear  scalar  transport.  Possibly,  the  most  important  features  of  Dyn-SGS  that  emerged 
from  this  study  are  the  following: 

•  For  smooth  problems,  this  model  does  not  deteriorate  the  nominal  order  of  accuracy  of  the 
spectral  approximation  of  the  governing  equations  (see  the  error  norms  computed  for  the 
geostrophically  balanced  flow  in  a  3D  channel.) 

•  It  is  flexible  and  robust  with  respect  to  the  flow  regime  and  grid  size. 

•  It  is  completely  free  of  a  user-tunable  parameter. 

•  In  the  neighborhood  of  sharp  gradients,  it  limits  and  controls  the  magnitude  of  the  over-  and 
under-shoots  without  compromising  the  solution  away  from  the  discontinuity.  Moreover,  the 
sharpness  of  the  discontinuity  is  very  well  preserved. 

It  was  also  shown  that  dynamic  stabilization  and  large  eddy  simulations  are  achieved  by  one  scheme 
alone. 

When  it  comes  to  parallel  performance,  the  cost  of  this  method  is  that  of  a  second-order  Laplace 
operator,  whose  computation  only  requires  one  communication,  against  the  two  (or  more)  necessary 
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Figure  14:  ID  Burgers  equation:  Characteristic  curves  of  the  computed  solution.  Left:  Ax  = 
6.25  x  10“3  m.  Right:  Ax  =  3.125  x  10“3  m. 
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Figure  15:  ID  Burgers  equation.  Point-wise  relative  error  in  logarithmic  scale  on  the  x  —  t  plane. 
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t  =  3/(271-) 


Figure  16:  ID  Burgers  equation  at  t  =  3/(2tt).  Global  solution  and  details  of  the  discontinuity 
(smaller  subplots).  Solution  computed  on  an  80  element  grid  of  order  8.  The  shock  is  well-resolved 
within  one  high-order  element.  The  over-  and  under-shoots  are  limited,  well  controlled,  and  do  not 
propagate  from  the  shock  across  O.  The  sharpness  of  the  discontinuity  is  well  preserved. 
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t  =  2/7T 


Figure  17:  Like  Fig.  16,  but  at  the  final  time  t  =  2/7T  s. 


x 


t  =  2  /  7T 


Figure  18:  ID  Burgers  equation.  Point-wise  error,  \u  —  uexact\,  at  the  two  time  levels  t  = 
[3/ (2tt),  2/7 r]  s. 
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Figure  19:  ID  Burgers  equation.  Structure  of  /z,  Eq.  (15),  at  the  three  time  levels  t  = 
[1/tt,  3/(2tt),  2/7 r]  s. 


Figure  20:  ID  Burgers  equation.  Normalized  L i,  L2,  and  Loo  error  norms  as  a  function  of  time. 
Left:  Ax  =  6.25  x  10“3  m.  Right:  Ax  =  3.125  x  10~3  m.  As  expected,  a  large  jump  in  the  error  for 
the  infinity  norm  occurs  as  soon  as  the  shock  forms. 


when  using  a  fourth  (or  higher)  order  hyper-diffusion.  On  the  way  towards  exascale  computing, 
fewer  communications  have  a  direct  impact  on  the  simulation  speed;  this  factor  is  of  fundamental 
importance  for  the  design  of  next  generation  atmospheric  models. 

We  have  seen  how  Dyn-SGS  is  based  on  element-wise  coefficients.  Based  on  the  observations  of 
[34;  3;  17],  the  current  method  (and  hence  the  solution)  would  benefit  if  the  discontinuous,  element¬ 
wise  viscosity  could  be  smoothed  via  some  proper  mechanism.  Taking  advantage  of  the  continuity 
of  the  solution  across  spectral  elements,  a  point-wise  definition  of  the  diffusion  coefficient  should  be 
sufficient.  This  issue  will  be  explored  in  the  future,  together  with  a  more  thorough  analysis  of  how 
Dyn-SGS  performs  on  passive  tracers  and  how  it  can  be  used  as  a  turbulence  model. 

We  have  anticipated  that  the  current  dissipation  is  independent  of  the  numerical  method.  It 
could  be  implemented  in  a  discontinuous  Galerkin,  finite  volume,  or  other  environments  alike.  As 
shown  in  the  Appendix,  we  have  been  working  on  its  application  using  discontinuous  Galerkin  and 
will  report  more  of  our  findings  in  a  future  paper. 
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A  Stabilized  discontinuous  Galerkin 

Dyn-SGS  is  independent  of  the  numerical  method.  Its  construction  is  only  tied  to  the  governing 
equations.  To  show  this,  we  applied  it  to  the  discontinuous  Galerkin  solution  of  the  rising  thermal 
bubble  problem  described  in  [20]  and  solved  therein  using  CG  and  DG.  These  results  are  reported 
in  this  Appendix  because  the  application  of  Dyn-SGS  to  the  discontinuous  Galerkin  method  falls 
beyond  the  scope  of  this  paper;  nevertheless,  it  is  meaningful  to  show  how  the  viscosity  described 
in  this  paper  can  be  utilized  outside  the  realm  of  finite  and  spectral  elements.  A  full  analysis  of  the 
performance  and  applicability  of  Dyn-SGS  to  DG  will  be  specifically  analyzed  in  a  future  work. 

For  DG,  the  flux  form  of  the  governing  equations  must  be  adopted,  so  that  we  re-write  equations 
(3)  in  conservation  form  as: 


dP  |  Pui  n 
dt  dxj 


dpuj  |  pUjUj  |  dp  _  drfjGS 
dt  dxj  dxi  dxj 

dj>8  pOuj  dQjGS 

dt  dxj  dxj  ’ 


pg$i, 


(34a) 

(34b) 
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Figure  21:  Rising  thermal  bubble  in  17  =  [1000] 2  m2.  DG  solution  at  A 9  at  i  =  700  for  the  two 
resolutions  Ax’  =  A z  =  12.5  m  (20  x  20  elements  of  order  4)  and  Ax  =  A z  =  6.25  m  (40  x  40  elements 
of  order  4). 


where  t^GS  and  QjGS  are,  again,  modeled  as  (4)  and  (6).  The  details  on  the  DG  approximation 
of  (34)  can  be  found  in  [33].  Here  we  simply  state  that  the  second-order  viscous  operators  are 
discretized  via  the  Local  Discontinuous  Galerkin  (LDG)  approach  [11].  To  build  Dyn-SGS  for  DG, 
the  equation  residuals  must  account  for  the  numerical  flux  that  results  from  the  DG  approximation. 
This  is  necessary  to  remove  the  inherent  dissipation  of  DG  from  the  SGS  model. 

The  solution  is  plotted  in  Fig.  21  for  two  grid  resolutions:  Ax  =  A^  =  12.5  m,  left,  and  Ax  = 
A z  =  6.25,  right,  in  the  domain  =  [1000] 2  m2.  In  [20]  the  stability  of  the  solution  was  achieved 
via  a  non-dissipative  low  pass  spatial  filter  [52;  6].  For  direct  comparison,  a  viscous  solution  of  the 
same  problem  is  reported  in  [53]  using  a  parameter-dependent  element-based  viscosity. 

We  are  fully  aware  of  the  many  other  stabilization  methods  designed  for  discontinuous  Galerkin 
(DG)  (see,  e.g.,  [44;  34]);  however,  the  fact  that  Dyn-SGS  is  based  on  the  physical  Navier-Stokes 
stress  operators  and  seems  feasible  for  large  eddy  simulations,  merits  consideration  for  proposing 
another  stabilizing  scheme  for  DG. 
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